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Abstract 

In this article, we present a state-of-the-art algorithm for solving the relativistic viscous hydro- 
dynamic equation with QCD equation of state. The numerical method is based on the second-order 
Godunov method and has less numerical dissipation, which are crucial in describing of quark-gluon 
plasma in high energy heavy- ion collisions. We apply the algorithm to several numerical test prob- 
lems such as sound wave propagation, shock tube and blast wave problems. In the sound wave 
propagation, the intrinsic numerical viscosity is measured and its explicit expression is shown, 
which is the second-order of spatial resolution both in the presence and absence of physical viscos- 
ity. The expression of the numerical viscosity can be used to determine the maximum cell size in 
order to accurately measure the effect of physical viscosity in the numerical simulation. 

PACS numbers: 
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I. INTRODUCTION 



A relativistic fluid approach is applied to various high energy phenomena in astrophysics, 
nuclear and hadron physics and has brought us a lot of interesting and outstanding results. In 
particular, recent relativistic hydrodynamic analyses revealed a new and interesting feature 
of Quark-Gluon Plasma (QGP) in high energy heavy-ion collisions. Since the Relativistic 
Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL) started its operation 
in 2000, a lot of discoveries have been made and gave us insight of quantum chromodynamics 
(QCD) phase transition and the QGP. One of the most interesting and surprising outcomes 
at RHIC was the production of the strongly interacting QGP (sQGP), which was confirmed 
by both experimental and theory sides. The highlights of them were (i) strong elliptic flow 
which suggests that collectivity and thermalization are achieved; (ii) strong jet quenching 
which confirms that hot and dense matter is created after collisions; (iii) the quark number 

n n 

scaling of elliptic flow which indicates that the hot quark soup is produced |2[. For 
these achievements, relativistic hydrodynamic models have made a lot of contribution. For 
example, at that time only hydrodynamic models could explain the strong elliptic flow at 
RHIC, which was considered as a direct evidence for production of sQGP at RHIC. Because 
of the success of relativistic hydrodynamic model at RHIC, at present hydrodynamic analysis 
is a useful and powerful tool for understanding dynamics of hot and dense matter in high 
energy heavy-ion collisions. 

In the early stage of the hydrodynamic studies at RHIC, viscosity effects were not taken 
into account. However, gradually, detailed analyses of experimental data in relativistic 
heavy-ion collisions reveal limitation of ideal hydrodynamic models. In Ref. jsj], for the first 
time, quantitative analyses of elliptic flow was done with relativistic viscous hydrodynamic 
model. They showed that ideal hydrodynamics overestimates the elliptic flow as a function 
of transverse momentum and hydrodynamic calculation with finite viscosity explains the 
experimental data better. Since then, main purpose of the phenomenological study for 
relativistic heavy-ion collisions at RHIC and LHC is sifted to obtain detailed information 
of bulk property of QGP such as its transport coefficients. Besides, recent high statistical 
experimental data at RHIC and LHC impose more rigorous numerical treatment on the 
hydro dynamical models. Recently both at RHIC and LHC the higher harmonic anisotropic 
flow which is the Fourier coefficient of particle yield as a function of azimuthal angle is 
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reported. One of important origins of the higher harmonics is considered as event-by-event 
fluctuations. To obtain the precise value of transport coefficients with relativistic viscous 
hydrodynamics, we need to choose an algorithm with small numerical dissipation and artifact 
for inviscid part carefully. Usually each algorithm has advantage or disadvantage in terms 
of coding, computational time, precision and stability. Up to now, unfortunately, only 
little attention has been paid to numerical aspects in hydrodynamic models for high energy 
heavy-ion collisions. 

In this article, we present a state-of-the-art algorithm for solving the relativistic viscous 
hydrodynamic equation with the QCD equation of state (EoS). Our physics application 
requires a numerical scheme that can treat a shock wave appropriately and has less numerical 
dissipation in order to gain comprehensive understanding of recent high energy heavy-ion 
collision physics. These advantages can be achieved by implementing a Riemann solver 
for the relativistic ideal hydrodynamics. In particular, we propose a new Riemann solver 
for the QCD EoS at low baryon density, which has not been considered in astrophysical 
application where baryon density is usually much higher. We derive our Riemann solver 
by analytically solving the relativistic Riemann problem in low baryon density within the 
approximation scheme proposed by 4| . As we will see in Sec. [V] where we perform several 
numerical tests, our new algorithm with the Riemann solver has an advantage over other 
prevailing algorithms such as KT, NT and SHASTA from the point of view of analyses 
for current relativistic heavy-ion collisions. By implementing our new Riemann solver for 
relativistic ideal hydrodynamics in a numerical scheme for causal viscous hydrodynamics 
recently proposed in Ref. ji]], we can also construct a new algorithm for causal viscous 
hydrodynamics for QGP. 

The article is organized as follows. In Sec. |TT1 we review current hydrodynamic models for 
relativistic heavy-ion collisions and introduce the basics of relativistic hydrodynamics. In 
Sec. IIIIl we explain the QCD EoS at high temperature and low baryon density based on the 
latest lattice QCD calculation. In Sec. IIVI we propose a new Riemann solver for the ideal 
fluid with the QCD EoS at high temperature and low baryon density. In Sec. IVJ using the 
numerical scheme, we show results of several numerical tests such as sound wave propagation, 
shock tube and blast wave problems. Section I VI I is devoted to summary and discussions. 
In this article, we adopt natural units; speed of light in vacuum c = 1, Boltzmann constant 
ks = 1 and Plank constant h = 1. 
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II. HYDRODYNAMIC MODELS 



First we list current hydrodynamic models which are applied to relativistic heavy-ion col- 
lisions 6| in Tables HI and HT1 Here we just mention important aspects of numerical simulations 
in relativistic hydrodynamic models. They are classified into ideal versions and viscous ones. 
One of important ingredients of hydrodynamic models is an EoS which is needed for solving 
the relativistic hydrodynamic equation and in which we can input various kind of physics re- 
lated to the QCD phase transition 1 . From comparison between hydrodynamic calculations 
and experimental data of high energy heavy-ion collisions, the information of the QCD phase 
diagram is obtained through the EoS used in the hydrodynamic calculation. The Bag model 
type EoS with the first order phase transition has been widely used in relativistic hydrody- 
namic models, because of its simplicity and lack of conclusive results on EoS of QCD. In 
recent hydrodynamical calculations, lattice(-inspired) EoS is gradually employed, because 
of the development of thermodynamical analyses based on the first principle calculation, 
lattice QCD simulation. 



TABLE I: Ideal hydrodynamical models. In the table, we use the following abbreviation. 1QCD: lat- 
tice QCD inspired EoS, SPH: smoothed particle hydrodynamics, PPM: piecewise parabolic method. 



Ref. 


dim 


EoS 


scheme 


Hama [7] 


3+1 


Bag model 


SPH 


Eirano [8] 


3+1 


Bag model 


PPM 


Nonaka [9] 


3+1 


Bag model 


Lagrange 


Hirano [10, 11] 


3+1 


1QCD 


PPM 


Petersen [12] 


3+1 


hadron gas 


SHASTA 


Holopainen [13] 


2+1 


1QCD 


SHASTA 



Another important ingredient in hydrodynamical models is numerical scheme for solving 
the relativistic ideal and viscous hydrodynamical equation. As long as we analyze multiplic- 



1 In addition to EoS, for the application to relativistic heavy-ion collisions further aspects should be dis- 
cussed; initial condition and termination condition (freezeout processes and final state interactions) of 
hydrodynamic simulation. Because these facts are more related physics in relativistic hydrodynamic 
collisions and out of scope of the manuscript, they are not mentioned in Tables HI and HT1 
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TABLE II: Viscous hydrodynamical models. In the table, we use the following abbreviation. CD: 
central difference, and KT: Kurganov-Tadmor (KT) scheme. 



Ref. 


dim 


EoS 


scheme 


Romatschke [3] 


2+1 


1QCD 


CD 


Luzum [14] 


2+1 


1QCD 


CD 


Schenke [15] 


3+1 


1QCD 


KT 


Song [16] 


2+1 


1QCD 


SHASTA 


Chaudhuri [17, .18] 


2+1 


Bag model 


SHASTA 



ities and collective flow using smoothed initial conditions, the choice of numerical scheme is 
not so important. However, when we start to investigate viscosity effects and event-by-event 
fluctuations in recent high statistic experimental data, we need to choose suitable numer- 
ical schemes carefully. For stability of hydrodynamic calculation, numerical dissipation is 
needed. Therefore in order to evaluate physical viscosity in high energy heavy-ion collisions, 
we need to avoid or control the artificial viscous effect in numerical calculation of relativistic 
viscous hydrodynamic calculation. One of candidates of accurate numerical schemes can be 
those with Riemann solvers for relativistic ideal hydrodynamics (references therein jlj]). The 
Riemann solver is a method to calculate numerical flux by using exact solution of the Rie- 
mann problems at the interfaces separating numerical grid cells, and can be used to describe 
the flows with strong shocks and sharp discontinuity stably and highly accurately. 

Here, we mention the basis of hydrodynamics briefly. The relativistic hydrodynamic 
equation is given by the conservation laws of energy, momentum, and baryon number: 

d fl T^(x) = 0, (1) 
= 0, (2) 

where T^ u (x) is the energy-momentum tensor and Jg{x) is the baryon current. In the case 
of relativistic ideal fluid, the energy-momentum tensor and baryon current are given by 

T^{x) = [e(x)+p(x)]u tl (x)u u (x) -p(x)g» v , (3) 
r B (x) = n B {x)u»{x) (4) 

where e(x), p(x), Ub(x), and u^(x) (u /i (x)w At (x) = 1) are the proper energy density, pressure 
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and baryon density which are evaluated in the rest frame of the fluid, and four velocity, 
respectively Eqs. (JT]) and fl2]) are solved numerically. 

When one starts to include the effects of dissipation into relativistic hydrodynamics, 
one is confronted with a rather complicated situation. One of the difficulties is that a 
naive introduction of viscosities, first order theory (i.e. first order in gradients) suffers from 
acausality. In order to avoid this problem, terms of second order in heat flow and viscosities 



have to be included in the expression for the entropy [19H25I]. but the systematic treatment 
of these second order terms has not been established. Although there is remarkable progress 
toward the construction of a fully consistent relativistic viscous hydro dynamical theory for 
description of relativistic heavy-ion collisions, there are still ongoing discussions about the 



26]. 



formulation of the equations of motion and about the numerical procedures 

At first order the new structures are proportional to gradients of the velocity field u M , 
and only three proportionality constants appear: the shear viscosity rj, the bulk viscosity 
(, and the baryon number conductivity a. At second order many more new parameters 
related to relaxation phenomena, such as relaxation times for each diffusive modes t v , t^, 
and T a appear. Currently, most of viscous hydrodynamical calculations use the relativistic 
dissipative equations of motion that were derived phenomenologically by Israel and Stewart 
[jjj which is chosen in this paper (see Appendix |A]) and variants of those 20 25| . Recently, 
a second-order viscous hydrodynamics from AdS/CFT correspondence was derived 271 ]. 
as well as a set of generalized Israel- Stewart equations from kinetic theory via Grad's 14- 



momentum expansion, which have several new terms [28||. On the other hand, however, a 
stable first-order relativistic dissipative hydrodynamical scheme was also proposed on the 
basis of renormalization-group consideration 



29 



30]. 



There are two choices for the local rest frame in relativistic viscous hydrodynamic equa- 
tion. One is Eckart frame 3l| where the direction of four velocity is the same as that of 



particle flux vector. The other is Landau frame 32[ where the direction of four velocity is 
the same as that of energy flux vector. Because in high energy collisions at RHIC and LHC 
baryon number density is very small (Sec. IIIII) . Landau frame is more suitable for QCD at 
high temperature and low baryon number density. 
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III. QCD EQUATION OF STATE AT HIGH TEMPERATURE AND LOW 
BARYON DENSITY 



In ultra-relativistic heavy-ion collisions at the LHC and RHIC, the relevant region in 
the QCD phase diagram is high temperature (T ~ 200-1000 MeV) and low baryon density 
(fiB ~ 0-100MeV) region. In this region, we expect that the QCD EoS can be approximated 
by taking into account leading order contribution of the finite baryon chemical potential. In 
other words, due to the charge conjugation (C) symmetry of QCD, the C-even quantities, 
e.g. pressure, energy density, temperature, and sound velocity, are approximated by those 
at vanishing baryon chemical potential, while the C-odd quantities, e.g. baryon density, and 
baryon chemical potential, are approximated by the first order contribution of the chemical 
potential. Note that in this approximation the C-even quantities are independent of /j,b, 
while the C-odd quantities depends on both T and \ib in principle. For example, 

p(T, (j, B ) = p(T, 0) + ± X (T, 0)/4 + 0(ji%) « p{T, 0), (5) 
n B (T, p B ) = MEl^A = x ( Tj o)^ + 0(ji%) « X (T, 0)/x fl , (6) 

where x stands for the baryon number susceptibility. Although the first principle lattice 
QCD simulation is limited at vanishing baryon chemical potential, we can access the ther- 
modynamic properties at low baryon chemical potential by using x(T, 0) in the above ap- 



proximation. Indeed, combining the result of the state-of-the-art lattice simulation [331. 134 1 
(p « 1-4 T 4 at T « 200-1000 MeV, x « 0.2-0.3 T 2 at T « 200-400 MeV) and the typical 
values of the baryon chemical potential in the heavy- ion collisions (fis ~ 24 MeV at RHIC 



and 1 MeV at LHC 



351]). we can estimate the importance of next-to- leading order term in 



the pressure by taking its ratio to the leading order term x{T,0)fi B /p(T,0) ~ 0.3(/xs/T) 2 , 
which yields only 0.4% correction at the RHIC and 0.00075% at the LHC. Therefore we 
regard this approximation is quantitatively successful in all the region of the QGP fireball 
at the RHIC and LHC. 

In numerical tests in Sec. |V] we will consider an EoS for free gas of gluons (free gas 
EoS) and that for realistic interacting quarks and gluons calculated by the lattice QCD 
simulation (lattice QCD EoS), which we plot in Fig. [TJ In the free gas EoS, we adopt the 
parameterization of 371 ]: 

48T 4 

e(T,0)=3p(T,0) = — , (7) 



7T 2 
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FIG. 1: (color online). Equation of states for free gas of gluons (free gas EoS) and interacting 
quarks and gluons calculated by the lattice QCD simulation (lattice QCD EoS) at vanishing baryon 
chemical potential. 

to make comparison with other numerical schemes and introduce x(T, 0) = eT 2 with e C 1 
in order to achieve effectively gluonic matter without quarks. In the lattice QCD EoS, we 
adopt the parameterization for the trace anomaly I = e — 3p for (2+1) flavors given in 
Eq. (3.1) and Table 2 of |33|: 



p(T,0) 



T 



dV J(T') 

rpl rp,4 



exp(-hi/t - h 2 /t 2 ) 



f ■ [tanh(/ 1 • t + h) + 1] 



(8) 
(9) 



1 + 9i ■ t + g 2 ■ t 2 

with t = T/(200 MeV), h = 0.1396, h x = -0.1800, h 2 = 0.0350, f = 2.76, /i = 6.79, 
f 2 = —5.29, 9\ = —0.47, and g 2 = 1.04. We parameterize the baryon number susceptibility 
X by fitting Fig. 7 and Table 1 of 34J: 

-T 



X(T,0) 
a 



aT 



1 + tanh 



AT 

0.15, T = 167 MeV, AT 



60 MeV. 



(10) 
(11) 



IV. RIEMANN SOLVER FOR THE IDEAL FLUID 
A. Exact solution of the relativistic Riemann problem 

Riemann problem is a classic one-dimensional initial value problem in hydrodynamics 
with infinitesimal dissipation and plays an essential role in numerical hydrodynamics. Since 
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we are interested in QCD matter in extremely high temperatures, we restrict our discus- 



sion to the relativistic hydrodynamics 



32]. The basic equations of the relativistic ideal 



hydrodynamics are the conservation equations for baryon number, momentum, and energy: 



d_ 

at 



hi 



\ E J 



( 



\ 



Dv 
mv + pi 
m 



\ 



J 



(12) 



where D, m, E are densities of baryon number, momentum, and energy; p, v are pressure 
and flow vector; and J is the identity matrix. The relation between the conservative variables 
U = (D, m, E) and primitive variables V = (ns, v,p) are 



D = 7n B , 

m = (e + p)~j 2 v, 



E 



(e + p)7 2 -P, 



(13) 
(14) 
(15) 



where 7 = (1 — 1 | 2 ) ™ 1//2 and e = e(p, ns) is given by the QCD EoS. 

The initial condition of the Riemann problem is given by two uniform states separated 
by a discontinuity surface at x = 0: 



V(x,t = 0) 



V L (x<0) 
V R (x > 0) 



(16) 



The exact solution to this problem is constructed from three types of self-similar flows, that 
is, flows which depend only on £ = x/t at t > 0: shock wave, rarefaction wave, and contact 
discontinuity 38] . Shock wave is a discontinuous surface moving at a constant velocity f s h, 
across which physical states are related by Rankine-Hugoniot jump conditions: 



[l/D] = -C [v x ] 

[(e + p)-fv x \ = C[p], 

[(e + jj)7%J = 0, 

[(e + p)7-p/D] = (W, 



(17) 
(18) 
(19) 
(20) 



where [q] = q — qs denotes the difference between the preshock state qs (S = L, R) and the 
postshock state q and ( = %h/j, 7sh = (1 — ^sh)" 1 ^ 3 = lshDs(v s h — v X; $)- Rarefaction wave 
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FIG. 2: An example of the exact solution for a Riemann problem. The regions L, L*, R* , and R 
are uniform and contact with each other on a shock wave, a contact discontinuity, and a rarefaction 
wave from the left. The region of the rarefaction wave is depicted by dotted lines. Whether L 
and L* (R and R*) are connected by a shock wave or a rarefaction wave depends on the relative 
pressure between L and L* (R and R*). 

is a continuous self-similar flow, through which physical states evolve by nonlinear ordinary 
differential equations: 

+ n B l Vy[v x - + ra B7 w 2 ( w x- - 0"^- = 0; ( 21 ) 

(e + p)7 2 K-0^ + (l-^)| = 0, (22) 
(e + p)7 2 K-0^f +%,e| = 0. (23) 

Contact discontinuity is also a discontinuous surface, across which pressure p and the flow 
velocity v x are continuous while other variables are discontinuous in general. 

The exact solution to the relativistic Riemann problem with general initial condition is 



given by Pons et al. 38( , which we summarize as follows. See also Fig. |2] as an example. 



1. Connect the initial left uniform state (L) and emergent uniform state inside (L*) by 
a shock wave or a rarefaction wave, depending on the relative pressure between these 
uniform states. 

2. Connect the initial right uniform state (R) and emergent uniform state inside (R*) by 
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a shock wave or a rarefaction wave, depending on the relative pressure between these 
uniform states. 

3. Connect the emergent uniform states (L*,R*) by a contact discontinuity. 

When the pressure of emergent uniform state is higher (lower) than the initial uniform state 
in contact, the shock (rarefaction) wave connects these uniform states. For the third step 
to hold, the emergent uniform states must be chosen so that the pressure p and the flow 
velocity v x is continuous across the contact discontinuity. 



B. Approximation scheme for the low-density QCD equation of state 



In numerical hydrodynamics with the shockwave capturing scheme, where each cell 
boundary is considered as independent Riemann problems, solving the ordinary differen- 
tial equation for the rarefaction wave Eqs. (l2ip -(l23 p numerically costs a lot of computational 



time. Mignone et al. (4] proposed an efficient approximation scheme to the exact solution 
of the Riemann problem. In their scheme, the rarefaction waves are approximated by the 
discontinuity that satisfies conservation laws. However, the original approximation scheme 
(4] needs to be modified for the QCD matter with low baryon density because it frequently 
uses the specific enthalpy h = (e+p)/n,B which diverges in vanishing baryon density tib = 0. 
Moreover, the model EoSs they have proposed, whose analytical simplicity also accelerates 
the numerical calculation, does not directly fit to the QCD EoS in low baryon density. 

Here we present a new approximation scheme for the QCD matter with low baryon 
density. First the Rankine-Hugoniot jump conditions Eqs. f [T7j) - (l2Q]) yield the Taub adiabat 



39] 



e +p x 
n B 



e(p; S) +p e s +p s 
n 2 B (p;S) n% s 



[P], (24) 



where V(p; S) (S = L, R) denotes the postshock variables, that is, the variables in the emer- 
gent uniform states (L*, R*). Once we specify a trial postshock pressure p and approximate 
the QCD EoS by e = e(p,n B ) ~ e(p,n B = 0) as given in Sec. HID 2 , we can solve the Taub 



2 If the baryon density is not low, one must solve e(p; S) = e(p, ub{p', S)) together with Eq. ([26]) to obtain 
e(p; S) and nsip] S), which requires additional iteration for solving this equation. 
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adiabat for the postshock variables V(p; S) with the preshock variables Vs, that is, the 
variables in the initial uniform states (L,R): 



e(p;S) = e(p,n fl = 0), (25) 

/ m / {e(p; S) +Ps}{e(p; S) + p} 
n B (p;S) = n B , s J {es + ps)(es+p) • (26) 



The baryon flux across the shock j(p; S) is also solved: 



[{e + p)/n%} 
2 e(p; S) + p s [p] 
HB ' S es+ps [e]-[ P Y (27) 



For later convenience, we define normalized baryon flux J(p; S) = j(p; S)/ub,s'- 

j2 fc5) = £feto W (28) 
es + Ps [ej - [pj 

with which the flow velocity w x (p; S) is given by 



(eg + Ps)ls v x,s + [p]C(p; 
{es + Ps)-fs + \p]{v x ,st(p) S) + l} 



v x (p;S) = - 2 ^ . r-.Tr . ,1 ^ , TT ' ( 2 9) 



C(p;5) = y r— 5 • (30) 

In practice, the limit [p] — >■ and [e] — > in the normalized baryon flux J(p; S) becomes 
numerically inaccurate. Therefore, when [p] or [e] is tiny, we switch to the analytical limiting 
value: 



b] ^o^'^ l- cl{p- S) l-c 2 sS 



l¥? JfoS) 2 = , ^Fzn =r Ig 4-. ( 31 ) 



where c s (p; 5) = c s (p,n B = 0) is the sound velocity. The sign in £(p; 5*) is chosen to be +(— ) 
for S = R(L) so that the correct shock propagation is ensured when the approximation 
scheme gives an exact solution. 

Since the two postshock states are separated by a contact discontinuity, pressure p and 
the flow velocity v x must be continuous. Therefore we have to solve 

v x (p;L) = v x (p;R), (32) 
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whose solution p* gives v x = v x (p*; L) = v x (p*; R) and other postshock variables. This part 
is solved numerically by the Newton-Raphson algorithm with following iteration: 

V (n + D = _(») _ v x & n hL)-v x (pM;R) 



v' x (p; S) 



dv x {p; S) 
dp 

{C(p; S) + [p]gp; S)}{1- v x ,sv x (p; S)} - v x ( P ; S) 
(e s + pshl + [p] K,sC(p; ^) + i} 



1 /'^W^,^ 



nfl ->Ul Aa^ " nB ' S 'd P V nj(p; 5) 



(34) 



M^°y ^MV-g)-v ' (35) 

d fe + p\ 2 d f e(p; S) + p 



GS + PS 'l- r eS ot P ^ (36) 



e(p; 5) + p 5 V e (p; 5) + p s (p; 5) 
This is the new approximation scheme for QCD matter with low baryon density. It is evident 
that there is no singularity in the limit % — > in the new approximation scheme. 



V. NUMERICAL TESTS 



By applying the Riemann solver in Sec. IIVI to the numerical scheme of causal viscous 
hydrodynamics 5], we solve several test problems, namely sound wave propagation, shock 
tube and blast wave problems in both ideal and viscous hydrodynamics. The structure of 
numerical algorithm is reviewed in Appendix [Bj 



A. Sound wave propagation 

Here we perform a simulation of sound wave propagation in ideal hydrodynamics using the 
numerical scheme presented above. The system length is L x = A = 2 fm and is discretized 
with iV ce ii = 50, 100, 200, and 400 cells. We set an initial condition 

V{x,t = 0)= (0, — -sin(27rx/A),0,0,po + ^sin(2vrx/A) J = V^x), (37) 

V c s0 (e +p ) / 

and impose a periodic boundary condition V(— A/2, t) = V(A/2, t). Here e = e(p , n B = 0), 
Cso = c s (po,riB = 0) and po = 10 3 fm -4 and bp = 10 _1 fm -4 . Since the amplitude of the 
wave is small 5p/po = 10 -4 <C 1, the nonlinear hydrodynamic equation is approximated by 
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FIG. 3: (color online). LI norm for pressure at t = X/c s q for the free gas EoS (left) and for the 



lattice QCD EoS (right). The dotted line indicates a scaling L(p;N cc \\) oc 1/N 



cell" 



linearized hydro-dynamic equation, which possesses a sound wave mode V s (x,t) = V; n i t (x — 
c s ot) as its solution. We analyze the precision of our numerical scheme and its dependence 
on TVceii by calculating the LI norm for pressure after one cycle t = X/c s q: 

L(p,p s ; TVceii) = ^2 I P(xi,\/c s0 ) -p s (xi,\/c s0 ) \ — — . (38) 



i=l 



N, 



cell 



We expect a scaling L(p,p s ; N ce n) oc (Sp/N^ cll ) ■ N ce]1 ■ (A/iV ccll ) = \5p/N^ cll after one cycle 
since our numerical scheme is of second order accuracy with respect to space and time 
discretization. The results of the LI norm for the free gas EoS and the lattice QCD EoS 
are shown in Fig. [3j We indeed find a scaling L(p,p s ; N cc \\) oc 1/A^ cll for both equations of 
states, which is consistent with the theoretical expectation. 

The simulation of sound wave propagation can also be utilized to estimate the numerical 
dissipation of the scheme. Since any numerical scheme introduces tiny numerical dissipation, 
the sound wave in the simulation is attenuated even without physical viscosity. The value 
of numerical dissipation is evaluated by the value of physical shear viscosity which gives the 
same amount of sound wave attenuation in the linearized region. By linear analysis, the 
disp ersion relation of the sound mode in viscous hydrodynamics with n B = 0, C = cr = 0is 



to 
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±c s0 k - i-fk 2 + 0(k 3 ), 
2r) 



(39) 
(40) 



3(e +Po)' 

Note that the dispersion relation is independent of the relaxation time for shear mode T n in 
long wavelength limit. The amplitude of sound wave with wave length A = 2n/k is decreased 
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FIG. 4: (color online). Numerical dissipation as a function of cell size Ax for (left) the free gas 
EoS and (right) the lattice QCD EoS. The dotted line indicates r/mim 

ps 1000(Ax) 2 for both EoSs. 



by a factor of exp 



8ir 2 ri 



3Ac s o(eo+po) 



after one cycle (t = A/c 



s0) 



p s (x, \/c s0 ; r})-po= \p s {x, A/c s0 ) - po] e 3Ac =o (e o +p o ) 



(41) 



To quantify the attenuation of the sound wave due to the viscosity, let us utilize the LI 
norm and define the numerical dissipation ?7 num : 

-A/2 



L(p s (ri),p s ;N cen oo) 



dx | p s (x, X/c s0 ; 7]) - p s (x, X/c s0 ) 



A/2 

2\Sp 



7T 



1 — e 3Ac so( e o+po) 



by which we obtain: 



3A 



c s o(e + p )ln 



7T 



2X5p 



L(p,p s ; iVcell) 



(42) 
(43) 



(44) 



In Fig. HI we show the numerical dissipation of our scheme for the free gas EoS and for 
the lattice QCD EoS. For both EoSs, the numerical dissipation can be approximated by 
Vmxm ~ 1000(Ax) 2 , where Ax = A/iV cc ii. From Eq. (14"4"1) and the second order accuracy of 
our numerical scheme L(p,p s ; N ce n) cx \5p/N^ eU = {Sp/X) ■ (Ax) 2 , the numerical dissipation 
is expected to scale with ?7 num oc [c s0 (e + Po)/X] ■ (Ax) 2 . Using the values of p , e , c s0 , and 
A, we find 



1 c s0 (eo + p ) (A ^ 2 
A 



(45) 



for both EoSs. 
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FIG. 5: (color online). (Left) LI norm L(p(r]),p s ; N cc \\) as a function of physical viscosity 
ij. The dotted line stands for the LI norm L\\ n (rj) in the linear analysis. (Right) LI norm 
L(p(rj),p s (r]); N ce n) as a function of physical viscosity rj. 

In Fig. [5j we show the numerical results of sound wave propagation for causal viscous 
hydrodynamics with the free gas EoS. In this calculation, we choose ( = a = and t v = 
10r]/sT and set an initial condition by Eq. ( 137|) . In the left panel, we show our numerical 
result of the sound wave attenuation due to both physical and numerical viscosities by 
calculating 

L(p(v)iPs;N ccU ) = ^ \p( x i^/ c so;v) - p s (xi,X/c s0 ) I — - . (46) 



i=l 



cell 



It shows that L(p(rj),p s ; N CQ \\) converges to ^1^(77) with larger iV cc u at fixed 77 and the con- 
vergence is faster at larger 77. This tendency is due to the numerical dissipation 77 num (Ax); 
when ?7 num (Ax) 77, the discretization effect is expected to be overwhelmed by the physical 
viscosity. In order to disentangle the physical and numerical viscosities, we calculate the 
following LI norm: 



L(p(v)iPs(v); N ccii) = ^2 I P( x *> V c so;?7) -Ps{xi, \/c s0 ;r)) 



A 



i=l 



N, 



cell 



(47) 



which eliminates the contribution from sound wave attenuation due to the physical viscosity. 
The result is shown in the right panel. We find that L(p(r)),p s (r)); N cea ) does not depend 
much on the physical viscosity 77. This indicates that the numerical dissipation 77 num w 
[c s o(eo +Po)/A] • (Ax) 2 gives a universal estimate of the numerical dissipation of our scheme 
in the presence of the physical viscosity. 3 



The increase of L(p(r]),p s (r]); N ce \i) with N cc \\ = 400 at large r\ > 2 fm 3 is due to the limitation of the 
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In our numerical scheme, the numerical dissipation is ?7 num m [c s (e + p)/\] ■ (Ax) 2 = 
(c s sT/X) ■ (Ax) 2 , where s denotes the entropy density. Since we are interested in physical 
viscosity of rj ~ (O.l-l)s, the condition (r) num /r)) ~ [c a T / (0.1-l)A]-(Ax) 2 <C 1 gives Ax <C 0.8- 
2.6 fm at T = 500 MeV and with A = 10 fm, which are the typical temperature and system 
length scale at the relativistic heavy-ion collisions. This condition becomes more severe at 
higher temperature or when finer structure is of interest. We emphasize that appropriate 
fine grid size calculation is indispensable for any physical observables in heavy-ion collisions, 
to discuss the value of physical viscosity from comparison with experimental data. 



B. Shock tube problem 

The shock tube problem is analytically solvable for perfect fluid with the free gas EoS. 
It provides an important test to know performance and accuracy of different numerical 
schemes. To compare our numerical algorithm to other numerical schemes (SHASTA, NT, 



KT schemes) and the analytical solution 36], we start the test calculation with the same 



initial conditions as those of Ref. [37]. The initial temperature on the left is To = 0.4 GeV 
and that on the right is To = 0.2 GeV. In the calculation we employ the free gas EoS. The 
spatial cell size and the Courant number are set to be Ax = 0.1 and A = 0.4 4 , respectively. 
Because a numerical calculation with enough fine grid and time step should converge on the 
analytical solution, the same discretization for spacial grid size and time step is important 
for accuracy test of numerical methods. 

Figure [6] shows the energy density distribution, the velocity and the invariant expansion 
rate with our algorithm, KT, NT and SHASTA together with the analytical solution. For 
these values, KT, NT, and SHASTA algorithms reproduce the analytical solution with almost 
the same accuracy and numerical artifacts. The difference between the analytical solution 



linear analysis in Eq. (|39[) . The dispersion relation has higher order contributions cj — ±c s ofc — ijk 2 ± 
{j/c s0 )(c 2 s0 t v - 7 /2)fc 3 + C(fc 4 ) [271. At 1 = V c s0, the third order term gives e ~ 0(1CT 4 ) • (r^/far 3 ) 2 



correction to the damped sound wave with Eq. (|39|) . which therefore is different from the exact solution 
by LI norm SL ~ 8p ■ e ■ A ~ 0(1O~ 5 ) • (?7/fm~ 3 ) 2 [fm~ 3 ]. Due to this difference, the LI norm Eq. (HT1) 
saturates at L > SL with large N ce i\. Similar argument holds for the limitation of the linear analysis for 
ideal hydrodynamics. In this case e ~ (Sp/po) ■ (l/c s o) ~ O(10 -4 ) and the LI norm Eq. ([55} saturates at 
L > SL ~ 5p ■ e ■ X ~ 0(1CT 5 ) [far 3 ] with large 2V ce ii. 
4 In our algorithm the Courant number is determined based on the Courant-Friedrichs-Lewvy (CFL) Con- 
dition, which realizes high-precise calculation (Appendix [B]) . 
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FIG. 6: (color online). The analytic (thine line) and numerical solutions of the relativistic Riemann 
problem on a grid with N x = 100 cells with Ax = 0.1 fm, after Nt = 100 time steps at t = 4 fm/c. 
(a) the energy density distribution e (b) the velocity v (c) the invariant expansion rate 9 with our 
algorithm (solid line) KT (dotted line), NT (dashed-dotted line) and SHASTA (dashed line). 



and numerical calculations indicates existence of artificial viscosity in numerical schemes. 
On the other hand, our numerical results are closer to analytical solution especially at 
x = 3 fm compared to KT, NT and SHASTA algorithms which suggests that our numerical 
algorithm contains less numerical dissipation. This tendency appears clearly in the invariant 
expansion rate 9 in Fig. [HI Only our numerical scheme follows the shape of the analytical 
solution from x = 3 fm to x = 5 fm. The numerical dissipation is indispensable for stability 
of numerical calculation of the relativistic hydrodynamical equation. However too much 
artificial viscosity smears numerical results and leads a wrong solution which is far from the 
analytical solution. 
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FIG. 7: Shear viscosity dependence of (a) energy density e, (b) velocity v and (c) invariant expan- 
sion rate 9. 

If an enough fine cell size is utilized in numerical calculation, the distinction among differ- 
ent algorithms becomes small, because numerical solution should converge to the analytical 
solution. However convergence speed to the analytical solution is different among numerical 
schemes. For example, to analyze the higher harmonics, numerical calculations with high 
statistics are needed, which indicates that we sometime reconcile a numerical calculation on 
coarser grids under current computational resources. According to physics target of appli- 
cation of hydrodynamics we need to choose the appropriate numerical method for solving 
the relativistic hydrodynamic equation. Besides, in relativistic heavy-ion collisions one of 
interesting and important topics is investigation of bulk property of the QGP such as its 
transport coefficients. To evaluate the physical viscosities of QGP from analyses of exper- 
imental data based on hydrodynamic models, we need to control the artificial viscosity in 
numerical calculation of hydrodynamical equation. The difficulty of distinguishing between 
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FIG. 8: (color online). EoS dependence of (a) energy density e, (b) velocity v and (c) invariant 
expansion rate 9. 

he physical viscosity and the artificial viscosity in numerical results was discussed in Ref. 



37[. For investigation of physical viscosity of QGP, the algorithm in which the artificial 
viscosity is well controlled is indispensable. 

In Fig. [7] shear viscosity dependence of energy density distribution, velocity and invariant 
expansion rate are shown. At finite shear viscosity, deviation from the analytical solution 
becomes large and shape of distribution is smeared. We observe the same tendency in finite 
bulk viscosity and baryon number conductivity calculation. 

Figure |8] indicates that EoS dependence of energy density distribution, velocity and in- 
variant expansion rate. For comparison, the same initial pressure distribution is employed 
for both cases. The fact that the sound velocity of lattice QCD EoS is smaller than that of 
free gas EoS (Fig. [T]) affects expansion rate. In Fig. [8] (c) the expansion rate of lattice QCD 
EoS is smaller that that of free gas EoS in almost entire region. As a consequence, velocity 
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of lattice QCD EoS is smaller than that of free gas EoS (Fig. [8] (b)) and expansion of shock 
wave in pressure distribution is smaller than that of free gas EoS (Fig. El (a)) . 

C. Blast wave problem 

We solve a (2+1) dimensional blast wave problem. The initial pressure and density are 
uniform and the initial flow vector is normalized to v r and points to the center of the system 



with pq = 1 fm -4 and v r = 0.9. The system area is a square with 6 fm x 6 fm and we 
discretize it with 384 points in each direction. 

We perform the blast wave simulation in (i) ideal and viscous hydrodynamics with free 
gas EoS and (ii) ideal and viscous hydrodynamics with lattice QCD EoS. In viscous hydro- 
dynamic simulations, we choose viscous coefficients r]/s = 0.1, ( = fm -3 , baryon number 
conductivity a = fm -1 , and relaxation time for shear mode t v = 10r]/sT = 1/T. 

In Fig. [9j we show the results of simulation (i) at t = 1.62 fm (1000 steps). In the upper 
panels, we plot the pressure and velocity profiles for the ideal hydrodynamic simulation. 
Note that the flow velocity field is dimensionless. At the center, we find a region with high 
pressure and vanishing flow velocity, which grows in time. In the lower panels, we show 
one- dimensional profiles of pressure and x(|/)-component of flow velocity at y(x) = fm for 
both ideal and viscous hydrodynamic simulations. It is clear that there is symmetry between 
x and y directions, which must be realized because of the initial condition. We find that the 
discontinuous change of pressure and flow velocity at \x\ ~ 0.7 fm in the ideal hydrodynamic 
simulation becomes continuous due to the finite shear viscosity. 

In Fig. [101 we show the results of simulation (ii) at t = 1.63 fm (1000 steps). In the upper 
panels, we plot the pressure and velocity profiles for the ideal hydrodynamic simulation and 
in the lower panels, we show one-dimensional profiles of pressure and x(?/)-component of 
flow velocity at y(x) = fm for both ideal and viscous hydrodynamic simulations. Here 
we find qualitatively same features as we find in the simulation (i), but we can find some 
quantitative differences. The pressure in the central region is about two times higher than 
that in the simulation (i). Also the size of the central region is about 10% smaller than that 
in the simulation (i). These differences are explained by the fact that the lattice QCD EoS 




(48) 
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FIG. 9: (color online). Simulation of ideal and viscous hydrodynamics with free gas EoS at 
t = 1.62 fm (1000 steps). The upper panels are two-dimensional profiles of (left) pressure and 
(right) flow velocity for the ideal hydrodynamic simulation. The lower panels are one-dimensional 
profiles of (left) pressure and (right) x(y)-component of flow velocity at y(x) = fm for both 
ideal and viscous hydrodynamic simulations. The finite viscous coefficient is r]/s = 0.1 and the 
relaxation time for the shear mode is = 1/T. 

is softer than the free gas EoS as shown in Fig. [TJ 

We have also successfully performed (3+1) dimensional blast wave simulations for both 



ideal and viscous hydrodynamics with the same initial condition p = 1 fm" 



0.9 as in 



the 2+1 dimensional simulations. In the viscous hydrodynamic simulation, we choose the 
same parameterization for viscosity t]/s = 0.1 and relaxation time t v — 1/T as before. Since 
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FIG. 10: (color online). Simulation of ideal and viscous hydrodynamics with lattice QCD EoS at 
t = 1 f m (550 steps). The upper panels are two-dimensional profiles of (left) pressure and (right) 
flow velocity for the ideal hydrodynamic simulation. The lower panels are one-dimensional profiles 
of (left) pressure and (right) x(y)-component of flow velocity at y(x) = fm for both ideal and 
viscous hydrodynamic simulations. The finite viscous coefficient is rj/s = 0.1 and the relaxation 
time for the shear mode is = 1/T. 



it is quite similar to the results of 2+1 dimensional case, we just remark it here but do not 
show the results. 
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VI. SUMMARY 



In this article, we presented a state-of-the-art numerical algorithm for solving the rela- 
tivistic viscous hydrodynamic equation with the QCD EoS. The numerical scheme is suitable 
for analyses of shock wave phenomena and has less numerical viscosity. Both features are 
important for understanding feature of QGP in high energy heavy-ion collisions. We apply 
the algorithm to several numerical test problems such as sound wave propagation, shock 
tube and blast wave problems. We investigated the precision of our numerical scheme in 
the sound wave propagation using the free gas EoS and the lattice QCD EoS. In both 
cases, LI norm scales as oc 1 / A^? cll with the number of cells, which shows the second or- 
der accuracy of our algorithm. Also, we estimated the numerical dissipation of our scheme 
Vmxm ~ [c a (e +p)/X] ■ (Ax) 2 both in the presence and absence of physical viscosity. We 
showed the results of shock tube test with our new numerical scheme which suffers less 
artificial dissipative effect. This suggests that it is more suitable for analyses of physical 
viscosities than SHASTA, KT, which are mainly used in studies of high energy heavy-ion 
collisions. We performed (2+1) and (3+1) dimensional blast wave simulations in ideal and 
viscous hydrodynamics with free gas EoS and lattice QCD EoS. 

The numerical scheme for relativistic viscous hydrodynamics with lattice QCD EoS is sta- 
ble, capable of capturing the shock wave and has less artificial viscosity. These features give 
us solid baseline for comprehensive understanding QGP in high energy heavy-ion collisions 
from the point of view of phenomenological analyses based on relativistic hydrodynamics. 
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Appendix A: Israel-Stewart equation 



Let us summarize the Israel- Stewart formalism for causal viscous hydrodynamics. In the 



32] 5 , energy- momentum tensor 



relativistic viscous hydrodynamics in Landau frame 
and baryon number current are decomposed as 

= eu»u v - (p + II) A"" + tt^, (Al) 
J% = n B u» + v%, (A2) 
N* v = f-uY, (A3) 

with bulk pressure II, shear stress tensor n^, and dissipative baryon number current u 1 ^ 
which satisfy ir^Uy = 0,ir^ = 0, v B u p = 0. Thermodynamic quantities e,p, and n B are re- 
lated through the equation of state p = p(e, n B ) derived in equilibrium state and the terms 
that include II, ir^, and u 1 ^ make extra contributions to T^ v and in non-equilibrium 
state. In the second order formalism by Israel and Stewart the entropy current s M is 
constructed so that it is defined locally without derivatives, includes terms with dissipa- 
tive quantities (II, z/^, ir^ u ) up to second order, and must satisfy the condition that s^u^ is 
maximized in equilibrium (II, v 1 ^, ir^ = 0). Such entropy current is then constructed by 

s» = su» - - ^(a n^ + am^u Bu ) - - fti/gi^ + f3 2 n^n pa ), (A4) 

with coupling coefficients «o,i and /3o,i,2 (> 0) and baryon number chemical potential fi B . 
In our algorithm, we neglect the couplings among different diffusive modes («o = ai = 0). 
Calculating the divergence of the entropy current by using the conservation laws Eqs. ([T|) 
and (El), we obtain 



J- * 7 

d^s" = - — (A^d^Uv + a d^ + O U) + —{d^u v - atd^u^ - fckfj) 

Tdfj, (^p) + aod^U + ai d w 7r; - (3 l u B ^ , (A5) 



where / = u^d^f. In deriving Eq. (1A5|) . we neglect terms in higher order deviation from 
equilibrium such as —^^-/3 Il 2 . In order to ensure that the entropy does not decrease, II, 



5 Since we are interested in the quark-gluon plasma at vanishing chemical potential, the Eckart frame, in 
which the four velocity u M is defined as Jg = ubu^ , is inconvenient. This is because the baryon number 
current is Jg = (0,0,0,0) in equilibrium and hence becomes ill-defined. 
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ir^, and Uq must obey the following constitutive equations: 



-n = c(A^v^ + «o<VB + A)ri), 



7T, 



2r]{{d^u u - axd^UBu - A^)), 



T / 



a d p Tl + ai9 CT 7r^ - ftz/, 



(A6) 
(A7) 
(A8) 



with bulk and shear viscosities (, 77 (> 0) and baryon number conductivity a (> 0). Here 
((X^)) denotes a spatial, symmetric, and traceless tensor extracted from a general tensor 

Xpcr ~\~ X a p 



A pcr A a/3 X Q/3 



(A9) 



2 3 

Note that the diffusive modes II, ir^, and are now dynamical degrees of freedom and 
relax toward the first-order Navier-Stokes values. The relaxation times for these diffusive 
modes are given by 



T"C = 0o(, t v = 2/3 2 ?7, T a = (3iG. 



(A10) 



Appendix B: Numerical algorithm 

Here we will give a brief summary of our numerical algorithm of the causal vis- 
cous hydrodynamics based on Ref. |5j. Let us first introduce the conserved variables 
U = (D, m, E) = C/;d + U v i S , where L/;d and U v i S denote contribution from the ideal and 
viscous components respectively. In the causal viscous hydrodynamics, we use the primitive 
variables V[d = (ub, v,p) for the ideal component and V V [ S = (n u , z/^) (z > j, i, j, k = 1, 2, 3) 
for the viscous component. Here we define as the total viscous component in the 
energy-momentum tensor IP" = —UA' 11 ' + ix^ v . Recall the relations among these variables 
[/id = C/id(^d) an d U vis = £/ vis (Vid, V v is)- The time evolution by the causal viscous hydro- 
dynamics is summarized by 

U(t) U(t + At), (Bl) 
V vis (t) -> V vis (t + At), (B2) 

where Eq. ( 1B1I) represents the time evolution by the continuity equations for the conserved 
quantities and Eq. (1B2I) represents that by the Israel- Stewart equation. In the numerical al- 
gorithm, we evolve the primitive variables V;d, v is according to these hydrodynamic equations. 
Note that U and "V^ is carry enough information to obtain U^, t/ vis , and V;d in principle. 
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Following Ref. jlj], we construct a numerical algorithm for the causal viscous hydrody- 
namics. We split the time evolution into 3 steps: 

1. Ideal part of the continuity equations: U id (t) ->■ U* d (t + At), V id (t) ->■ V£(t + At). 
Here we evolve only Ui d (t) by using the currents of the ideal component. The Riemann 
solver we proposed is used in evaluating the numerical flux to obtain the currents. In 
the ideal hydrodynamics, we only need this step. In the viscous hydrodynamics, this 
step does not give Ui d and V;d at time t + At as indicated by asterisks *. Note 
that in this step the total conserved quantities vary from U(t) = U[ d (t) + U v [ S (t) to 
U* d (t + At) + Uvis(t) so that this step satisfies the conservation law. 

2. Israel- Stewart equation: V vis (t) ->> V vis (t + At/2) ->> V vis (t + At). 

This step consists of solving the convection part and the relaxation part of the Israel- 
Stewart equation separately and implementing them by the operator splitting method 



as in Ref. 



a 



In the relaxation part, we calculate the time differentiation of Vi d in the 
Navier-Stokes terms by [V£(t + At) - V id (t)} /At. We use V£(t + At) to evaluate V id 
in all the other parts in the Israel- Stewart equation. We use the intermediate state 
Vvi s (t + At/2) in the next step. 

3. Viscous part of the continuity equations: U* d (t + At) + U vis (t) — > U{t + At). 

Here we evolve the sum U* d (t + At) + U v i S (t) by the currents of the viscous component 
using V vis (t + At/2) and V£(t + At) to obtain U(t + At). From U(t + At) and 
Vyi s (t + At), we get U id (t + At) ^ U* d (t + At) and V id (t + At) ^ V^(t + At) in general. 
Note that this step also satisfies the conservation law. 

The At is defined so that it satisfies the CFL condition of the telegrapher equation as 
in Ref. When the spatial dimension is more than 1, all the steps but the relaxation 
part in 2 are implemented with the dimensional splitting method. This numerical algorithm 
is applicable to both Landau and Eckart frames of causal viscous hydrodynamics. By this 
numerical algorithm, we achieve the second order accuracy in both space and time as we 
checked in Sec. [V] 
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FIG. 11: (color online). Temperature of the shock tube problem at t = 2.7 fm (1400 steps) in 
Landau and Eckart frames. The particle diffusion in Landau frame, or the heat conductivity in 
Eckart frame, is k = 0.05 fm -2 . 

Appendix C: Comparison of Landau and Eckart frames 



Here we compare simulations of viscous hydrodynamics in Landau and Eckart frames. 
The purpose of this comparison is to show that the original hydrodynamic code in Eckart 
frame |5[ is extended correctly to a code in Landau frame. Therefore we perform the simu- 
lation with the equation of state with high baryon density as in jfj] and we do not use the 
Riemann solver that we propose in the text. 

We simulate a shock tube problem with an initial condition: 



V(x,t = 0) 



V L = {duo, 0, v y0 , 0, p L0 ) (x < 0) 
V R = (dm, 0, -v y0 , 0,p R0 ) (x > 0) 



(CI) 



with dio = Plo = 10 fm -4 , dRo = pro = 1 fm -4 , and v y0 = 0.2. Note that d(x) denotes mass 
density in this simulation. The system length is 2 fm and is discretized with 200 cells. The 
particle diffusion in Landau frame, or equivalently the heat conductivity in Eckart frame, is 
k = 0.05 fm -2 and shear and bulk viscous coefficients are rj = £ = fm -3 . 

Shown in Fig. [TT] is the temperature at t — 2.7 fm (1400 steps) in Landau and Eckart 
frames. According to 19[, the difference in thermodynamic quantities in these frames are 
small. Since we cannot find frame dependence in the temperature profiles, we conclude that 
we have successfully extended the original code to the one in Landau frame. 
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